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Nonlinear  Least  Squares 
Paraaeter  Optimization 


A es  re  ACT  fCuM—  a»  rurmm  «aa  W — 1>  1 ease  Ay  me  —Up 

Marquardt's  algoritha  is  an  iaproved  technique  for  the  least  squares  estima- 
tion of  nonlinear  paraaeters.  The  present  report  describes  the  algorithm  and 
documents  a FORTRAN  subroutine  for  carrying  it  out.  A specific  application 
of  the  subroutine  is  discussad  in  datsil;  fitting  e complex  yew  equation  to 
data  obtained  fron  a Free  Flight  Range.  An  example  is  given  in  which  the 
fitting  process  converged  to  the  wrong  answer  without  the  algorithm  and  to 
the  right  answer  with  it. 
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I . INTRODUCTION 

Fitting  an  equation  to  & boss  of  data  points  (x,  y)  is  an  ancient 
and  honorable  pastime.  Circa  1800,  Gauss  and  Legendre  discovered  inde- 
pendently what  proved  to  be  a giant  advancement  in  technique:  the 

Method  of  Least  Squares.  Refinements  of  this  basic  method  are  still 
showing  up  in  the  literature.  This  report  concerns  one  such  refinement, 
proposed  by  D.  W .Marquardt 1 . 

Suppose  re  are  given  an  equation  of  the  form 

y - f U.  P)  (1) 


where  P is  a set  of  n unknown  constant  parameters: 


We  are  also  given  a set  of  m measurements 


Y 


(2) 


(3) 


at  the  corresponding  (presumably  error-free)  x values  xjt  x^,  . . . x^, 

where  m exceeds  n:  there  are  more  measurements  than  unknown  parameters. 

For  convenience,  wc  assume  each  measurement  y^  is  equally  reliable;  this 

eliminates  the  need  for  weighting  factors,  which  add  nothing  to  the 
concepts  involved. 

Wc  are  asked  to  determine  the  least  squares  fit  of  Equation  (1) 
to  the  set  of  measurements.  That  is,  we  seek  those  values  of  the  a^'s 

that  minimize  c,  the  sum  of  the  squares  of  the  residuals  of  the  fit: 

m 

e * F (P)  2 l [y  - f (x.,  P)]2  (4) 

i-1  1 1 


The  criterion  function  F in  (4)  is  assumed  to  be  an  analytic  func- 
tion of  the  n parameters.  We  can  describe  the  situation  geometrically 
as  follows.  Consider  an  n-diroensional  Euclidean  space  S in  which  any 
set  of  parameter  values  constitutes  the  coordinates  of  a point  P.  Then 


1 . Donald  W.  Marquardt,  "An  Algorithm  for  Least-Squares  Estimation  of 
Nonlinear  Parameters,"  J.  Sop,  Indust.  Appl.  Math.,  V ol.  11,  No.  2, 
pp  431-441,  June  1363. 
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Preceding  page  blank 


e is  the  value  of  the  continuous  scalar  point  function  P at  point  P.  Por 
each  point  P in  the  parameter  space  S the  .re  corresponds  a single  value  of 
e and  our  task  is  to  search  S in  sons  well-defined,  automated  Banner  for 
the  point  (and  hopefully  there  is  only  one  such  point)  that  yields  the 
absolute  mini bub  value. 

In  practice,  the  criterion  function  often  possesses  one  cr  store 
local  minimum  values  in  addition  to  the  desired  absolute  minimum.  The 
particular  minim*  we  reach  then  depends  on  where  in  spac<  3 we  start 
our  search.  If  we  start  too  close  to  the  point  L associated  with  oue 
of  these  local  ainiaa,  and  if  our  aatheaatical  searching  procedure  — in 
the  interest  :vf  computer  econo ay  — is  not  prohibitively  sophisticated, 
then  we  ere  going  to  be  drawn  irresistibly  to  L.  Even  if  we  start  far 
froa  L,  a funny  thing  can  happen  on  the  way  to  the  right  answer:  our 
path  nay  pass  through  L*s  sphere  of  influence  and  we  are  trapped.  We 
will  give  an  exaaple  of  this  misfortune  in  the  final  section  of  this 
report. 

A necessary  (but  insufficient)  condition  for  e to  be  at  a local  or 
absolute  niniaua  is  that  the  gradient  of  c be  the  zero  vector: 


Thus  we  seek  the  set  of  paraaeter  values  that  satisfy  all  n component 
equations  of  (S)  simultaneous ly « 

When  Equation  (1)  is  linear  in  the  parameters: 

y • ^ U)  ♦ *2  *2  (x)  + . . . + an  *n  (x)  (6) 

Equation  (5)  represents  a system  also  linear  in  the  parameters . In 
principle,  such  a system  (provided  it  is  nonsingular  and  strble)  can  bu 
solved  by  any  one  of  a variety  of  well-known  techniques.  This  is  too 
easy;  thus  we  specify  here  that  Equation  (1)  is  nonlinear  in  at  least 
one  of  the  parameters. 

For  this  nonlinear  situation,  two  of  the  estat  ished  methods  for 
minimizing  e are: 

(a)  steepest  descent  (alias  gradient  search); 

(b)  differential  corrections  (alias  Taylor  series  expansion, 
alias  Gauss -Newton  method) . 

Maxquardt  has  proposed  a third  technique.  As  we  will  see,  his 
algorithm1  represents  a blend  of  the  first  two  methods,  retaining  the 
best  features  of  each. 
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In  Sections  II  to  IV,  we  will  discuss  enough  of  the  differential 
corrections  end  steepest  descent  techniques  to  see  what  is  involved  in 
■n  interpolation  of  the  two.  In  Sections  V to  VII,  we  will  discuss  the 
Marquardt  algorithm  and  its  iaplementation  as  a FORTRAN  subroutine  MARQ. 
In  Sections  VIII  to  XI,  we  will  show  how  subroutine  MARQ  is  used  to  fit 
a complex  yaw  equation  to  a set  of  measurements  obtained  from  on  enclosed 
'>>©» flight  spark -photography  range.  In  the  final  section  we  will  com- 
pare a cose  run  with  and  without  the  old  of  Marquardt' s algorithm. 


II.  DIFFERENTIAL  CORRECTIONS 

Ne  con  obtain  the  basic  equations  of  the  differential  corrections 
technique  by  approximating  y or  c or  grad  e by  a first-order  Taylor 
expansion  about  a given  point  PQ.  For  our  purposes,  the  simplest  choice 

is  grad  e: 

*-N.*  £ [“N."’ 


where 


Aa^  «*  the  change  in  a^  in  moving  from  PQ  to  P 

and  where  the  zero  subscript  denotes  evaluation  at  point  PQ.  Under 

approximation  (7) , the  nonlinear  system  (5)  is  replaced  by  a system  that 
is  linear  in  the  parameter  increments  Aa^ : 


M*  "[^1 


, k « 1,  2,  , . . n (8) 


or,  simplifying  the  notation, 


t ... 


Aa.  ■ e,  , k ■ 1,  2. 


.Ml' 


(11) 


and  where  the  factor  1/2  has  been  introduced  to  simplify  later  expres- 
sions for  0^  and  otjk. 

The  system  (9)  can  be  written  in  matrix  form  as 

7T  • a (12) 


where 

a - the  n x n symmetric  matrix 

AP  • the  row*  vector  |Aa  , Aa  , . . . Aa„  1 

V i 2 njb 

£ * the  row*  vector  ^0^,  8^,  . . . 8nj^ 

The  symmetric  matrix  a is  sometimes  called  the  curvature  matrix  because 
it  is  a measure  of  the  curvature  of  the  c hypersurface  in  the  parameter 
space.  From  (10),  we  see  that  I is  a vector  in  the  direction,  of  the 
negative  gradient  of  e at  Pq: 


? “ - y («**d  e)( 


(13) 


Equation  (9),  or  the  equivalent  Equation  (12),  represents  a system 
of  n linear  equations  in  the  n increments  Aa^ , Hence  we  can  solve  the 

system  for  these  increments: 


Aa. 


k-1 


6k  “jk*  J ■ !»  2*  • • • n 


(14a) 


or 


A*  -1  • * 


(14b) 


*If  we  had  written  the  product  in  Eq.  (12)  as  a • AP,  then  AP  and  £ 

would  be  oolum  veotore. 


! 

i 

1 

! 


> 


j 
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where 


“■(^4 


the  inverse  of  matrix  a 


ajk 


cofactor  of  element  a 
determinant  of  a 


2Jl 


Of  course,  to  solve  the  system  (9),  we  must  be  able  to  express  8^ 
and  in  terms  of  the  given  fitting  function  of  Hq.  (1).  Substituting 
Kq . (4)  in  (10),  we  have: 


-T« 


- 2 


A faf  (x.,  P)1 

Z Cyi  " f CV  p0)]  5rH 

i-1  L Jc 


m 

Z - f <V  poJ] 


ik 


(IS) 


i-1 


where 


ik 


(xA,  P)] 


(16) 


Similarly,  substituting  (4)  in  (11),  we  have 
m 


ajk  * 


Z { ■>«  “lk  - - f <v  p0)]  [-4y  } »7> 


The  second-order  term  in  (17)  could  be  a big  nuisance  if  we  let  it,  but 
traditionally  — and  with  ample  justification  — • this  term  is  dropped. 
Hence  our  working  definition  of  is  the  approximation 


“jk  “ 


z ^ 

i-1 


(18) 
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When  the  system  (12)  has  been  solved  for  the  increment  vector,  we 
can  then  make  the  differential  corrections  on  PQ»  that  is,  we  can  then 

obtain  a new  point  whose  vector  is 


P„  ♦ LP 
0 


(19) 


and  whose  associated  error  value  is  e ■ F (P  ) . 

1 1 

If  Eq.  (7)  were  exact  (which  would  happen  only  if  the  original 
fitting  function  in  Eq.  (1)  were  linear  in  the  parameters  a^),  then  P 

would  be  the  desired  point  at  which  e is  a minimum.  Since  Eq.  (7)  is 
only  an  approximation,  P is  not  the  solution  to  the  nonlinear  problem. 

Indeed,  if  the  starting  point  PQ  is  not  close  enough  to  the  solution 

point  (and  it  is  usually  difficult  to  say  beforehand  how  close  is  "close 
enough") , then  e will  be  larger  than  e0;  the  process  may  or  may  not 

recover  from  this  inauspicious  start.  By  the  rules  of  the  game,  P is 

made  the  new  starting  point,  whether  or  not  it  is  an  improvement. 

System  (12)  is  then  re-solved  for  a new  set  of  parameter  increments, 
which  in  turn  yield  — for  better  or  for  worse  — » a new  point  P^  and 

so  on,  each  iteration  of  the  scheme  taking  U3  to  the  next  point.  This 

process  may  converge  to  a unique  point  PM  and  if  so,  = F (PM)  may  be 

the  absolute  minimum  value  of  the  criterion  function.  But  don't  count 
on  it. 


An  estimate  of  the  error  of  the  fit  at  any  point  P is  usually 
obtained  by  the  relation 


E *»  b 


(20) 


where 


b * 1.0  to  obtain  an  estimate  of  the  root -mean -square  error 
(more  precisely,  the  RMS  deviation  between  the  measure- 
ments and  a hypothetical  set  of  "true"  measurements 
free  from  observational  errors) 

» 0.67449  to  obtain  an  estimate  of  the  probable  error 
and  an  estimate  of  the  error  in  any  parameter  a.  at  point  P is  given  by 

J 

Ej  ■ E (21) 
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Because  of  its  application  in  Equation  (21),  the  inverse  matrix  a is 
sometimes  called  the  error  matrix. 

Before  we  consider  the  steepest  descent  approach,  it  will  be  con 
venient  to  simplify  whatever  physical  dimensions  are  involved  in  our 
equations . 


III.  SCALING 


If  the  parameters  are  dimensional  and  — worse  yet  — not  all 

of  the  same  dimensions,  then  our  parameter  space  S can  be  a hodge-podge 
of  femto-drams,  fermis  and  fortnights.  Letting  [ denote  "dimensions 

of",  we  have 


• K]d 


1 


Md  ■ [y2/a;  ‘i 

H.  ■ N.  ■ [yH 


(22) 


Suppose  now  that  we  introduce  a scaled  parameter  space  S in  which 


(23) 


so  that 


[y]d»  k * 1,  2,  . . . n 


(24) 


Then  the  correspondingly  scaled  partial  derivatives  and  a,  $ and  grad  e 
elements  are,  by  Equations  (16),  (18)  and  (15): 


11 


with  dimensions 


Moreover,  Cauchy's  inequality  applied  to  definition  (18)  asserts  that 

“jk2  4 “jj  v 


so  that 


- 1 < aj]{  < 1 (27) 

Thus  we  have  transformed  into  a space  § in  which: 

(a)  each  parameter  and  each  component  of  the  vectors  grad  e and 
has  the  same  dimension,  namely  [y]^; 

(b)  "ci  is  a dimensionless  matrix  whose  diagonal  elements  are  all 
unity  and  whose  off-diagonal  elements  satisfy  Inequality  (27) . 

All  equations  of  the  previous  section  hold  for  the  transformed 
elements;  in  particular.  Equation  (9)  becomes 

n 

£ ajk  • AS.  * k = 1,  2,  . . . n (28) 

j = l 

with  solution 

n 

= X)  \ ^jk’  j ■ 2*  • • • n (29) 

k=l 

Transformation  (23)  is  often  used  in  linear  least-squares  to  im- 
prove the  numerical  behavior  of  the  computations:  in  particular,  in- 

version of  the  curvature  matrix.  For  our  nonlinear,  iterative  problem, 
the  transformation  could  have  the  same  salutary  effect  but  at  the  cost 
of  a little  more  work.  This  is  because  the  values  of  the  scale  factors 
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depend  on  the  current  set  of  parameter  values.  Hence  each  tine 
the  parameters  are  up-dated,  a new  transformation  must  be  made:  a new 

A 

S space  created.  This  is  not  a big  problem. 

The  transformation  to  ^ is  not  essential  in  the  differential  cor- 
rections method;  it  is  needed  more  in  the  steepest  descent  approach. 
Some  of  the  properties  of  the  latter  technique  are  not  scale-invariant, 
so  that  choosing  the  right  space  can  be  important.  At  the  very  least, 
it  is  convenient  to  proceed  by  steepest  descent  in  a space  in  which  the 
components  of  grad  e are  dimensionally  equal. 


IV.  STEEPEST  DESCENT 

Provided  that  |grad  e|  is  not  zero  at  the  current  point  PQ  (if  it 

were,  PQ  would  be  the  desired  solution),  then  -grad  e at  that  point 

is  a vector  in  whose  direction  t will  decrease  most  rapidly  (at  least 
at  first)  as  we  move  away  from  PQ.  If  Pg  is  any  other  point  along  this 

negative  gradient,  we  have  from  Eq.  (13): 

T (h)  *T0  * ht  (30) 

where  h is  a dimensionless  positive  constant  and  where  the  three  vectors 
are  to  be  resolved  into  components  in  the  scaled  space  s . Thus  in 
component  form,  (30)  becomes: 


- h (31) 

In  the  steepest  descent  technique,  we  choose  that  h for  which  e is 
a local  minimum  along  the  gradient.  Perhaps  the  simplest  way  to  do  this 
is  to  sample  c at  appropriately  small  intervals  as  we  move  away  from  PQ 

along  the  negative  gradient.  As  soon  as  we  reach  a point  where  e has 
increased,  we  interpolate  between  that  point  and  the  previous  sampling 
point  to  determine  where  c has  a local  minimum.  Then  we  evaluate  the 
negative  gradient  at  this  new  point  and  start  off  again  in  the  new 
direction. 

The  difficulty  with  this  approach  is  that  in  the  neighborhood  of 
the  solution  point,  where  |grad  e|  is  nearly  zero,  further  progress  is 
painfully  slow.  Often,  the  sampling  size  must  be  shortened  beyond  all 
endurance.  Ingenious  variations  on  the  basic  steepest  descent  theme 
have  lessened  this  difficulty  but  not  removed  it. 
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V.  INTERPOLATION 


Comparing  the  two  methods  just  discussed,  we  note  that: 

(a)  Far  from  the  solution  point,  the  steepest  descent  technique 
is  superior!  It  must  proceed  so  as  to  decrease  c,  whereas  the  differ- 
ential corrections  method  has  no  such  obligation  and  is  likely  to 
diverge. 

(b)  Close  to  the  solution  point,  the  differential  corrections 
method  is  superior.  It  converges  in  the  very  region  where  the  steepest 
descent  technique  languishes . 

Marquardt1  has  proposed  an  interpolation  between  the  two  methods: 
a technique  that  behaves  like  the  steepest  descent  when  we  are  far  from 
the  solution  and  like  the  differential  corrections  method  when  we  enter 
the  neighborhood  in  which  the  linear  truncation,  Equation  (7) , is  ade- 
quate . 

To  achieve  this  interpolation,  a positive,  non-dimensional  constant 
X is  added  to  each  diasonal  element  or  the  transformed  curvature  matrix. 


at  is.  Equation  (28)  is  replaced  by 


0)  As  X ->■  0,  Equation  (32)  approaches  Equation  (28),  that  is. 
Equation  (32)  reduces  to  the  differential  corrections  net hod. 

The  interpolation  sends  us  from  a given  point  PQ  in  the  direction 
of  the  vector 


satisfying  (32),  whereas  the  gradient  method  sends  us  in  the  direction 
of  the  vector  The  angle  between  these  two  directions: 


always. lies  between  0 and  90*.  (This  angle  is  the  sane,  of  course,  in 
S and  5.)  Maxquardt  proves  that  for  the  given  point  PQ  (and  hence  for 

a given  I),  6 is  a continuous,  monotonically  decreasing  function  of  X, 
such  that: 

(a)  As  X -»  0 -*■  0. 

(b)  As  X -»  0,  0 approaches  some  value  9^  less  than  90°. 

(Marquardt  states  that  6^^  usually  lies  between  80*  and  90*,  indicating 

that  the  differential  corrections  nethod  usually  proceeds  almost  at 
right  angles  to  the  gradient  approach.  However,  in  the  yaw  fit  appli- 
cation to  be  discussed  later,  1 have  found  that  8^  is  usually  less 

than  80*  and  sometiaes  less  than  30*.) 

The  insertion  of  X into  (28)  makes  interpolation  possible;  an 
effective  interpolation  is  achieved  only  when  we  have  a sound  strategy 
for  changing  the  value  of  X according  to  whether  e has  increased  or 
decreased  since  its  previous  evaluation.  In  the  next  section,  we 
present  the  strategy  proposed  by  Marquardt. 


VI.  MARQUARDT 'S  ALGORITHM 

Let  e « F (P  ) be  the  value  of  the  criterion  function  at  some 

W O ^ 

current  point  PQ.  For  any  X,  we  can  solve  (32)  for  the  increments  AP, 

obtain  P - P,  + AP  and  evaluate 
1 0 

- F (Pj)  - G (P0.  X) 
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where  the  second  functional  fora  emphasizes  the  dependence  of  on 
both  the  previous  point  PQ  and  the  current  value  of  X. 

In  the  Marquardt  approach,  the  point  is  only  a candidate  for 
the  point  that  will  replace  PQ  as  our  current  set  of  parameter  values. 

A new  rule  is  this:  we  won't  move  from  PQ  to  any  new  point  P unless  e 

is  smaller  at  the  new  point  than  at  PQ: 


For  any  point  PQ  that  is  not  already  a local  minimum,  there  always 

exists  (at  least  in  theory)  a XM  such  that  for  X > XM,  (36)  holds.  The 

question  is:  how  much  larger  than  XM  should  we  choose  our  X?  Marquardt 's 

answer:  no  larger  than  necessary  within  a factor  of  A > 1.  That  is,  if 

our  initial  X doesn't  work,  try  AX,  A2X,  A3X,  etc.,  quitting  as  soon  as 
(36)  holds.  Ten  is  the  suggested  value  for  A. 

Of  course,  we  could  determine  X so  as  to  minimize  e locally,  just 
as  with  h in  the  steepest  descent  approach.  Marquardt  points  out,  how- 
ever, that  this  would  usually  result  in  a much  larger  X than  needed  to 
satisfy  (36);  we  would  be  faced  with  all  the  disadvantages  of  the  steep- 
est descent. 

On  the  other  hand,  we  would  like  to  simulate  the  differential 
corrections  method  as  soon  as  we  get  within  the  radius  of  convergence 
of  that  technique.  Hence,  when  things  are  going  well,  we  want  to  reduce 
X,  say  by  factors  of  A again.  Marquardt  suggests,  however,  that  if  X is 
already  negligible  compared  with  1 to  the  number  of  significant  digits 
carried,  then  there  is  no  need  to  reduce  X further. 

Marquardt 's  strategy  can  then  be  sunmarized  by  the  flow-chart  of 
Figure  1.  (This  is,  at  any  rate,  my  understanding  of  Marquardt 's 
strategy;  his  presentation*  differs  considerably  in  manner  and  form.) 

Note  from  the  figure  that  in  practice  Marquardt  replaces  (36)  by 
the  condition 


< eQ  (36a) 

Hie  distinction  may  seem  academic  but  it's  not.  In  an  early  version  of 
our  program  for  implementing  the  algorithm,  I unwisely  used  (36)  instead 
of  (36a).  Why,  after  all,  should  we  move  to  P^  if  the  error  there  is 

I very  soon  found  out  why  and  the  answer  (upon  after- 


the  same  as  at  P ? 

o 


• Ao  * A v 4 

j using  and  Absolve  E002)  for  ap 

• COMPUTE  P?*  A"P 

• COMPUTE  €,  • F (P,  ) 


thought)  is  obvious.  It  is  possible  for  the  current  point  PQ  at  the 

end  of  some  iteration  to  "be"  the  solution  point  in  the  sense  that  the 
machine  will  be  unable  to  find  a point  with  a smaller  associated  error 
(to  the  number  of  digits  carried).  In  this  situation,  of  course,  we 
should  exit  the  loop  triumphantly  and  proceed  tc  analyze  results.  The 
trouble  is:  any  convergence  test  for  determining  whether  or  not  to 

start  another  iteration  is  necessarily  imperfect.  The  test  is  usually 
based  solely  on  how  much  the  current  point  PQ  and/or  error  eQ  differs 

from  the  values  of  the  previous  iteration.  Thus  we  could  close  in  on 
the  answer  so  quickly  that  the  solution  point  differs  from  the  previous 
point  by  more  than  our  criterion  allows.  In  this  event,  another 
iteration  will  begin:  a hopeless  search  for  a candidate  P satisfying 

(36) . 

As  an  example  of  what  can  happen,  consider  the  test  case  whose 
results  are  summarized  in  Table  1.  The  fitting  equation  here  was  the 
complex  yaw  equation,  to  be  discussed  later:  the  convergence  criterion 
was  that  (1  - e^/e Q ) be  less  than  0.00001.  In  the  first  iteration,  the 

process  obtained  an  e (.005307)  less  than  eQ  (.007418),  so  that  it 

would  make  no  difference  whether  condition  (36)  or  (36a)  was  used  in 
the  scheme  of  Figure  1.  The  convergence  condition  was  not  satisfied 
(since  1 - t^/c Q * 0.285)  so  the  second  iteration  was  began.  Similarly, 

the  third  and  fourth  iterations  were  begun.  It  was  only  in  attempting 
to  complete  the  fourth  iteration  that  the  superiority  of  (36a)  to  (36) 
because  evident. 

The  condition  (36a),  < cQ,  was  finally  met  (in  this  case,  when 

achieved  equality  with  eQ  to  the  sixteen  decimal  digits  carried*)  at 

X = 1013.  When  condition  (36),  < eQ,  was  substituted  for  (36a),  the 

fourth  iteration  was  never  completed;  program  execution  was  interrupted 
at  X = 10 16  by  a floating  point  overflow.  The  hang-up  probably  occurred 
in  evaluating  the  determinant  of  matrix  y,  Gq.  (33).  When  X is  large, 

k 

say  X = 10  where  k > 2,  then  the  determinant  is  approximately  the  prod- 
uct of  the  diagonal  elements  of  y: 


DET  a (1  + X)n  a 10kn 

If  a computer  can't  handle  real  numbers  larger  than  10’1,  then  we  are  in 
trouble  when  k exceeds  j/n.  For  our  computer,  j - 154  and  for  the  test 
case,  n ■ 10;  hence  the  hang-up  at  k ■ 16. 


*0f  oouree,  it  is  unlikely  that  either  eQ  or  is  correct  to  the  six- 
teenth digit ; the  vagaries  of  round-off  error  alone  may  make  the  last 
feu  digits  meaningleee . 
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Test  Case  Using  the  Marquardt  Algorithm  of  Figure  1 with  A 
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In  our  United  experience,  the  test  case  of  Table  1 is  not  typical ; 

we  actually  contrived  this  unimpressive  exanple.  Still,  it  could  have 

happened  in  "real  life",  so  let  us  try  to  profit  by  it. 

We  see  from  Table  1 that  a lot  of  needless  work  was  done  in  the 

fourth  iteration  — the  matrix  equation  (32)  was  solved  twenty  times  — 
only  to  finish  at  the  point  already  obtained  in  the  third  iteration. 
Clearly,  there  is  room  here  for  improvement  in  technique. 

Of  the  various  ways  of  reducing  the  number  of  calculations,  the 
simplest  might  to  be  to  introduce  two  acceptable  conditions: 

e < e (36) 

l o 

and 

ef)  < < cQ  (1  + A)  (36b) 

where  A is  assigned  some  small  positive  value.  If  condition  (36)  is 
satisfied,  we  up-date  the  current  point  and  its  associated  error  and 
proceed  as  in  Figure  1.  If,  instead,  condition  (36b)  is  satisfied,  we 
assume  that  the  point  reached  at  the  conclusion  of  the  previous  iteration 
is  as  good  as  we  are  going  to  get.  h this  event,  the  entire  process 
ends:  that  previous  point  and  its  associated  error  are  taken  to  be  the 

final  point  and  error.  For  example,  if  we  take  A « 0.00001,  condition 
(36b)  would  be  satisfied  at  the  end  of  the  eighth  loop  — iteration 
4(8)  — of  Table  1.  At  this  stage,  the  program  would  assume  that  the 
solution  point  is  the  point  obtained  by  the  third  iteration. 

The  above  technique  may  be  an  excellent  practical  way  out  of  our 
difficulties;  yet  there  is  something  vaguely  unpleasant  about  condition 
(36b).  In  effect,  we  are  giving  up.  The  suspicion  would  always  linger 
that  in  another  loop  or  two  within  the  iteration,  the  situation  might 
suddenly  improve  and  condition  (36)  would  be  satisfied. 

We  ask  then:  how  can  we  cut  down  on  the  computations  without  relax- 

ing our  standards;  that  is,  how  can  we  satisfy  (36a)  with  less  effort? 

One  obvious  answer  is  suggested  by  Eq.  (34): 

AP  - £/X,  X » 1 (34a) 

a relationship  confirmed  by  Table  1.  For  large  values  of  X,  it  is  no 
longer  necessary  to  obtain  new  AP  candidates  by  the  time-consuming 
process  of  increasing  X and  re-solving  Eq.  (32)  by  matrix  inversion. 

All  succeeding  AP  candidates  within  that  iteration  will  be  very  nearly 
parallel  and  hence  can  be  obtained  directly  by  a scale  change: 

- ^rejected  /*>  * > 1 <37> 

candidate  candidate/ 
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where  S need  not  have  the  ease  value  assigned  to  A. 

We  now  ask:  Is  it  necessary  to  wait  until  X is  large  to  replace 
(32)  by  (37)?  Marquardt  Motions  that  in  some  circumstances  (charac- 
terized by  matrix  elements  exceeding  0.99)  he  found  it  helpful  to 

revert  to  Eq.  (37)  as  soon  as  the  angle  6 fell  below,  say  45*  (that  is, 
when  he  was  within  45*  of  the  steepest  descent  direction).  Usually  open 
to  a good  suggestion,  I incorporated  this  little  sub-scheme  in  the  grand 
plan  as  a feature  operative  under  all  circiastances . That  is,  whenever 
6 is  less  than  some  specified  value  G and  condition  (3£a)  isn't  satisfied 
for  vector  P ♦ &P,  the  next  vector  considered  is  P ♦ AP/B. 

The  final  fora  of  our  modified  Marquardt  algorithm  is  shown  in 
Figure  2.  If  the  user  feels  that  the  9 < G feature  is  too  great  a 
violation  of  the  spirit  of  Marquardt 's  algorithm,  he  need  only  set  G 
at  zero  to  avoid  that  feature. 

Table  2 shows  the  results  of  applying  the  modified  algorithm,  with 
A ■ B ■ 10,  G ■ 45*,  to  the  test  case  of  Table  1.  The  first  three 
iterations  are  identical,  but  the  fourth  has  been  improved  considerably. 
The  reduction  in  the  number  of  loops  in  the  fourth  iteration  is  not  as 
important  as  the  fact  that  each  loop  required  significantly  fewer 
machine  calculations.  Of  course,  we  could  introduce  condition  (36b)  into 
this  new  scheme  as  well.  With  A ■ 0.00001,  the  process  would  then  end 
at  iteration  4(2) . 

After  we  have  discussed  our  yaw  equation,  we  will  present  another 
test  case  in  more  detail.  First,  however,  we  introduce  the  FORTRAN  sub- 
routine MARQ  for  carrying  out  the  algorithm  of  Figure  2. 


VII . COMMENTS  ON  SUBROUTINE  MARQ 

At  least  five  programs  to  implement  Marquardt' s algorithm  pre-date 
the  FORTRAN  subroutine  MARQ  listed  in  Appendix  A: 

(a)  A FORTRAN  program  "Least -Squares  Estimation  of  Nonlinear 
Parameters"  cited  in  Reference  1 and  available  as  IBM  Share  Program 
No.  1426. 

(b)  A program  "Non-Linear  Least  Squares  (GAUSHAUS),"  University 
of  Wisconsin  Computing  Center,  written  by  D.  A.  Meeter  in  1964. 

(c)  A "Nonlinear  Least-Squares  Curve  Fitting  Program,"  a 1966 
modification  of  (b)  above  by  F.  S.  Wood.  A User's  Manual  but  no  list- 
ing of  this  program  is  given  in  Reference  2;  the  program,  a FORTRAN 


2.  Cuthbert  Daniel  and  Fred  S.  Wood,  Fitting  Equations  to  Data:  Computer 
Analysis  of  Multifaotor  Data  for  Soientistsand  Engineered  Wiley-  ~ 
Intern aienoe,  New  Io*kt  1971. 
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Table  2.  Same  Test  Case  as  in  Table  1*  but  Using  the  Modified  Marquardt  Algorithm 
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listing,  the  User's  Manual  and  test  problems  are  available  from 

(1)  SHARE  Library 
COSMIC 

University  of  Georgia 

Athens , Georgia  30601 

(Ask  for  Program  No.  360D-13.6.007) 

and  from 

(2)  VIM  Library 

Software  Distribution  Department 

Control  Data  Corporation 

3145  Porter  Drive 

Palo  Alto,  California  94304 

(Ask  for  Program  No.  G2-CAL-NLWOOD) 

(d)  A FORTRAN  subroutine  CURFIT  by  P.  R.  Bevington,  listed  and 
discussed  in  Reference  3. 

(e)  A FORTRAN  program  of  the  same  title  as  (a)  above,  cited  in 
Reference  3 and  available  as  IBM  Share  Program  EID-NLIN  No.  3094.01. 

An  early  version  of  our  MARQ  subroutine  was  based  on  Bevington's 
CURFIT3  --  to  which  I am  indebted  --  but  even  our  first,  tentative 
efforts  differed  from  CURFIT  in  various  minor  ways.  (Among  certain 
programmers  myself  included  — the  temptation  to  tinker  with  a 
presented  program  is  well-nigh  irresistible.  The  Law  of  Mutual 
Superiority4  applies:  "Anything  you  have  programmed,  I can  program 

better;  anything  I have  programmed,  you  can  program  better.")  Our 
current  MARQ  is  barely  recognizable  as  a descendant  of  CURFIT. 

Subroutine  MARQ  has  12  arguments: 

(EQS,  X,  Y,  M,  N,  P,  C,  E,  D,  R,  TH,  EK) 

Here  and  throughout  this  paper  the  usual  convention  applies  to  FORTRAN 
real  and  integer  variable  names:  integer  names  and  only  integer  names 

start  with  I,  J,  K,  L,  M or  N.  The  first  five  arguments  of  MARQ  are 
inputs : 


3.  Philip  R . Bevington , Data  Reduction  and  Error  Analuaia  for  the 
Physical  Sciences , McGraw-Hill,  Inc. > New  York,  1969. 

4.  Richard  V.  Andree,  Josephine  P.  Andree  and  David  D.  Andree,  Computer 

Proqramting : Techniques,  Analysis  and  Mathematical  Prentice-Hail 

Series  "in  Automatic  Computation , Hew  Jersey , 1973. 
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EQS  * the  dummy  name  of  the  subroutine  that  computes , for  a given 
parameter  point  P,  the  arguments  E,  D and  R defined  below. 
The  actual  name  used  when  calling  MARQ  must  be  declared  in 
on  EXTERNAL  statement  in  the  program  that  calls  MARQ.  (In 
our  yaw  analysis  example,  the  actual  nasie  used  is  EYAN. 
Details  of  subroutine  EQS/EYAW  will  be  discussed  in  Section 
X.) 

X,  Y ■ data  vectors  (xp,  (y^)  where  y^  is  the  observed  value  of 
the  dependent  variable  y at  x ■ x^. 

M ■ the  number  of  data  points  (x^,  yp 
■ the  dimension  of  vectors  X and  Y 

N » the  number  of  parameters.  Note:  if  N exceeds  10,  the 

dimensions  of  arrays  ALPHA,  BETA,  GAMMA,  PB  and  S on  line 
MARQ  3 (see  Appendix  A)  must  be  increased  to  equal  or 
exceed  N and  the  argument  10  on  line  MARQ  48  must  be 
increased  to  equal  the  number  of  rows  declared  for  GAMA. 

The  next  five  arguments  (P,  C,  B,  D and  R)  serve  as  both  input  and  out- 
put; the  last  two  arguments  (TH  and  BK)  are  solely  outputs.  Before  we 
discuss  these  arguments  individually,  it  will  be  helpful  to  see  the 
over-all  pattern  os  indicated  in  Table  3: 


Table  3.  Contents  of  MARQ  Arguments  P,  C,  E,  D,  R,  TH  and  EK 
when  Marquardt's  lambt’a  is  Used 
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From  this  table  we  note  that  of  the  last  seven  arguments,  only  P and 
C Dust  contain  inputs  the  first  tine  MARQ  is  called.  This  first  call 
merely  obtains  the  error  E and  arrays  D and  R associated  with  the  initial 


point  P . Thereafter,  each  call  returns  a new,  improved  point;  the  values 


returned  in  P,  C,  E,  D and  R after  the  K-th  call  serve  as  input  for  the 
(K  + l)-th  call.  We  now  discuss  the  arguments  in  detail: 


P 3 the  current  point  P 


that  is,  the  current  set  of  parameter 

In  each  appli- 


values  { a_j ) in  whatever  units  are  convenient. 

cation,  before  MARQ  is  called  the  first  time,  array  P must 
contain  the  initial  estimated  values  of  the  N parameters  in 
the  chosen  units.  This  same  set  is  returned  from  the  first 


call;  on  subsequent  calls,  the  input  P is  replaced  by  an 
improved  P upon  return. 


C = 


a flag  (initially)  and  Marquardt's  X thereafter.  In  each 
application,  before  MARQ  is  called  the  first  time,  C must 
be  set  to  - 2.0  or  - 1.0.  (The  fact  that  C is  negative 
alerts  MARQ  that  it  is  being  called  for  the  first  time.) 
Caution:  use  a FORTRAN  name,  not  the  number  - 2.0  or  - 1.0, 

in  the  argument  list,  since  MARQ  changes  the  value  of  this 
argument . 


Set  C = - 2.0  if  the  intent  is  to  ignore  the  X factor 
and  use  a straight-forward  differential  corrections  technique. 
(This  option  is  helpful,  for  example,  if  the  user  wants  to 
compare  the  process  with  and  without  Marquardt's  X for  a given 
set  of  data.)  Upon  return  from  the  first  call  of  MARQ,  C will 
have  the  value  zero  (representing  a zero  X)  and  will  remain  at 
zero  through  subsequent  iterations. 


Set  C = - 1.0  if  the  intent  is  to  use  the  Marquardt 
algorithm.  Upon  return  from  the  first  call  of  MARQ,  C will 
have  the  value  0.001A  (representing  A times  the  initial 
value  of  X).  For  the  second  and  subsequent  MARQ  calls,  the 
output  C will  be  the  value  of  Marquardt's  X that  was  required 
to  obtain  the  improved  point  concurrently  stored  in  array  P. 
We  have 


Output  C = A • (Input  C) 


where  r is  some  integer  equal  to  or  greater  than  - 1 (so  that 
the  smallest  possible  — and  usually  the  most  felicitous  -- 
output  C is  (l/A)-th  the  input  C) . 


E = 


a measure  of  the  error  of  the  fit  at  the  point  concurrently 
stored  in  array  P.  This  measure  is  computed  in  the  subroutine 
EQS  called  by  MARQ  and  hence  the  coder  of  EQS  defines  this 
argument,  preferably  as  the  RMS  or  probable  error  given  by 
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Eq.  (20).  MARQ  will  always  return  an  E value  less  than  or 
equal  to  the  input  E. 


the  matrix  of  partial  derivatives,  Eq.  (16),  evaluated  at 
the  point  concurrently  stored  in  array  P: 

D(I,  K)  - D.v 


3ak  x - X (I) 

I P ■ current  point 

a the  vector  of  residuals  of  the  fit  at  the  point  concurrently 
stored  in  array  P: 


R (I)  - Y (I)  - f (x,  P) 


x - X (I) 

P ■ current  point 


TH  » the  angle  8,  Eq.  (35),  in  degrees.  Upon  return  from  the  first 
call,  6 is  zero.  Thereafter,  8 is  the  angle  between  the  nega- 
tive gradient  at  the  input  point  (not  necessarily  the  best  way 
to  go)  and  the  straight  line  from  input  point  to  output  point. 
This  angle  was  available  within  MARQ  and  simple  curiosity 
prompted  me  to  make  it  an  output  argument. 

EK  a the  vector  of  estimated  errors  in  each  of  the  N parameters 

stored  in  array  P,  in  the  same  units  as  the  parameters.  These 
estimates  are  related  to  output  argument  E by  a modified  form 
of  Eq.  (21)  and  hence  if  E is  the  RMS  or  probable  error  of  the 
fit,  then  array  EK  contains  the  RMS  or  probable  error  estimates 
for  the  parameters. 

In  the  MARQ  listing.  Appendix  A,  note  that  the  three  parameters  A, 

B and  G of  Figure  2 are  defined  in  a DATA  statement,  line  MARQ  4 (where 
they  are  named  AN,  BN  and  GN,  respectively).  Thus  the  user  can  easily 
change  the  arbitrary  values  shown,  subject  to  the  conditions 


A > 1, 


B > 1 


and  the  suggested  limitation 


0 < G < 45  degrees 

The  feature  of  using  Eq.  (37)  when  9 ' G can  be  eliminated  by  defining 
GN  to  be  zero  in  the  DATA  statement. 
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The  data  in  Tables  1 and  2 were  obtained  with  the  aid  of  some  ad 
hoc  WRITE  statements  within  MARQ;  these  statements  do  not  appear  in  the 
list' ng.  Appendix  A. 

In  addition  to  subroutines  MARQ  and  EQS,  a matrix  inversion  sub- 
routine (called  on  line  MARQ  48)  must  be  included  in  the  over-all  pro- 
gram. We  have  used  an  available  subroutine  MATINV  listed  in  Appendix  B. 
(The  user  may , of  course,  substitute  the  inversion  scheme  of  his  choice, 
modifying  the  CALL  MATINV  statement  as  necessary.)  Upon  return  from 
MATINV,  the  matrix  GAMMA  is  replaced  by  its  inverse.  MATINV  contains 
one  output  statement : a SINGULAR  MATRIX  reprimand,  which  is  printed 

only  when  something  is  very,  very  wrong.  The  integer  6 in  the  state- 
ment WRITE  (6,  17)  may  have  to  be  changed  to  specify  the  proper  output 
unit  at  the  user's  installation. 

Subroutines  MARQ,  MATINV,  EYAW  (our  actual  name  for  the  dummy  EQS) 
and  EY  (our  program  that  calls  MARQ)  were  written  for  use  on  ARDC's 
BRLESC  I and  BRLESC  II  computers5.  Although  the  compilers  in  these  two 
computers  implement  a FORTRAN  that  is  not  quite  "standard"  and  not  quite 
FORTRAN  IV  (and  indeed  not  quite  the  same  on  the  two  computers) , it  is 
believed  that  the  above-mentioned  subroutines  have  been  restricted  for 
the  most  part  to  standard  statements  and  features.  One  exception  is 
the  nonstandard  function  ARCGOS  used  in  MARQ,  line  62,  to  determine  0. 
The  user  can  replace  ARCGOS  (TR)  with 

ATAN  (SQRT  (TR  **  (-  2)  - 1.0)) 


if  necessary. 

An  important  omission  in  the  programs  should  be  noted:  no  provision 

was  made  for  double-precision  arithmetic.  In  BRLESC  I/II,  all  real 
numbers  are  carried  with  a precision  of  sixteen  decimal  digits  (which  is 
double  precision  on  many  computers)  and  the  DOUBLE  PRECISION  type  decla- 
ration, though  permitted,  acts  only  as  a REAL  declaration.  For  computers 
where  double  precision  isn't  the  norm,  the  following  DOUBLE  PRECISION 
type  declarations  are  recommended: 

in  MARQ  : GAMMA 

in  MATINV  : A,  T1 

in  EYAW  : S 


6.  Lloyd  V.  Campbell  and  Glenn  A.  Beokt  BRLESC  I/II  FORTRAN,  Aberdeen 
Reeearoh  and  Development  Center  Teohnioal  Report  No.  5,  AD  704343 , 
Mardh,  1970. 
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VIII.  THE  YAK  EQUATION 

For  alssiles  fired  in  either  of  the  two  Free  Flight  Spark  Photography 
Ranges6*7  operated  by  the  Exterior  Ballistics  Laboratory,  we  can  obtain 
at  each  spark  station  a measurement  of 

(a)  the  down-range  distance  z,  metres 

and  (b)  the  dimensionless  complex  yaw  { » + i 

(For  our  purposes  — to  illustrate  the  use  of  Marquardt's  algorithm  — 
definitions  of  the  coordinate  systems  involved  are  of  no  interest.) 

Basically,  the  problem  is  to  fit  a given  complex  yaw  equation 

? - f (z,  P)  (38) 


to  a set  of  measurements 

{(V  ?Hi*  (z2*  ?H2*  5V2D*  * * ‘ (zm*  W 5Vm)} 

where  the  z measurements  are  considered  error-free  (though  of  course 
they  aren't).  For  most  (say  99.44%)  of  the  rounds  fired  in  the  two 
ranges,  Equation  (38)  is  assumed  to  have  the  form 


C - K e1*!  + K o1*2  + £D 
1 2 R 


(39) 


where 


Kj  - Kjo  e\>(z  " Zo) 

- Aj  + Bj  (z  - zQ)  + Cj  (z  - zQ)2 


}- 


1,  2 


(40) 

(41) 


SR  - yaw  of  repose 
g (B  + B ) 

3 - i iL 


B B V 2 
1 2 0 


(42) 


6.  Valter  K.  Rogers,  Jr.,  'The  Transonic  Free  Flight  Range,"  Ballietio 
Res  ear oh  Laboratories  Report  No.  1044,  June  1958,  AD  200177 . 

7.  Valter  F.  Braun , " The  Free  Flight  Aerodynamics  Range , " Ballietio 
Research  laboratories  Report  No.  1048 , July  1958,  AD  202249. 
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The  known  constants  are 


zQ  ■ reference  z (usually  aid-range) , m 


V ■ velocity  at  z . 
o o 

g « acceleration  of 

■ 9.80  a/sec2  for 

and  the  ten  unknown  parameters 

P(l)  - n ■ in  K 
1 10 

P(2)  - Xx,  1/a 
P(3)  « rad 
P(4)  ■ B , rad/ a 
P(5)  ■ C , rad/m2 


m/sec 

gravity 
t]he  two  ranges 

are 


P(6) 

■ n * In  K 

2 : 

P(7) 

- X 1/m 
2 

PCS) 

■ A , rad 
2 

P(9) 

■ B , rad/m 
2 

P(10) 

- C , rad/m' 
2 

(43) 


It  might  seem  more  straightforward  to  define  P(l)  and  P(6)  to  be  the 

arm  lengths  K and  K , respectively.  However,  pre-Marquardt  experi- 
10  20 

ence  has  indicated  that  for  our  yaw  problem,  the  differential  corrections 
process  is  a little  more  likely  to  diverge  when  the  K^o's  are  handled 

directly  as  parameters.  For  one  thing,  AKjQ  may  be  so  poor  that  (unless 

constraints  are  added),  KjQ  can  go  negative.  This  difficulty  doesn't 

arise  with  Hj : any  value  of  yields  a positive  K^.  Of  course,  in 

practice  our  starting  point  PQ  is  usually  so  well -determined  that  it 

would  make  no  difference  whether  we  worked  with  the  rj ^ * s or  the  K^'s . 


Complex  expressions  and  complex  arithmetic  are  not  allowed  in  the 
BRLESC  I/I I FORTRAN;  hence  we  must  separate  Eq.  (39)  into  its  real  and 
imaginary  parts: 


eH  ■ fH  (-z>  p)  S KJ  cos  + K2  cos  + 5r 

* fV  (Z»  P)  3 \ Sin  + K2  Si" 


(44) 


obtaining  two  fitting  equations . 
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IX.  COMMENTS  ON  HANDLING  TWO  FITTING  EQUATIONS 


Previously  in  this  report  we  have  been  considering  a single  fitting 
equation.  Eq.  (1).  To  handle  two  fitting  equations,  as  in  (44),  we 
proceed  as  follows.  Suppose  that  NL  is  the  nuaber  of  spark  stations  at 
which  measurements  were  taken.  Then  the  value  of  the  error  function,  Eq 
(4),  is  obtained  by  the  relation 

NL 

- I 

i«l 


{>«Hl 


- fH  (*i( 


p)]2  + [5 vi  " £V  (zi» 


} 


(45) 


Since  MARQ  has  no  provision  for  inputting  more  than  one  dependent 
variable,  we  obtain  the  equivalent  of  Eq.  (45)  by 

(a)  setting  M,  the  number  of  measurements,  at  twice  NL, 


and  (b)  defining 


X (I)  « z4  - z. 


Y (I)  - £ 


Hi 


Y (I  + NL)  ■ £, 


Vi 


1,  2,  . . NL 


There  is  no  need  to  define  the  second  half  of  vector  X's  M components, 
since  it  would  only  duplicate  the  first  half.  (Note  that  we  need  a 
minimum  of  six  spark  stations  in  order  for  M,  che  number  of  measurements 
to  exceed  ten,  the  number  of  parameters.)  Equation  (45)  then  becomes 

NL 

e - l [R  (I)2  + R (I  ♦ NL)2]  (46) 

1-1 


where  the  M residuals  are  obtained  by  the  relations 

R (I)  = Y (I)  - f (X  (I),  P)  1 

" P I ■ 1,  2,  . . NL  (47) 

R (I  + NL)  « Y (I  + NL)  - fy  (X  (I),  P)  J 

The  only  remaining  change  required  to  handle  two  fitting  equations 
instead  of  one,  occurs  in  the  formation  of  the  partial  derivative 
matrix  D.  Since  there  are  two  equations  and  ten  parameters,  twenty 
partial  derivatives  must  be  defined,  ten  of  the  form 

8 fH  (X  (I),  P) 




D (I,  K) 


(48) 


and  ten  of  the  fora 


3 fv  (X  (I),  P) 

D(I  ♦ NL,  K)  - Vrk5 (49) 

To  summarize,  the  changes  needed  to  handle  two  fitting  equations  occur 
in: 


(a)  the  manner  of  defining  M,  X and  Y in  the  program  that  calls 

MARQ 

and  (b)  the  manner  of  defining  E,  R and  D in  the  subroutine  EQS  dis- 
cussed in  the  next  section. 

The  technique  for  dealing  with  two  fitting  equations  can  easily  be 
extended  to  any  number  of  simultaneous  fitting  equations.  Note  that  in 
our  example,  and  are  dimensionally  equal.  When  the  dependent 

variables  are  not  all  of  the  same  dimensions,  a little  additional  work 
is  needed  to  insure  dimensional  consistency  in  Eq.  (46)  and  elsewhere. 


X.  COMMENTS  ON  SUBROUTINE  EQS/EYAW 

The  user  must  code  the  subroutine  whose  dummy  name  is  EQS.  The 
arguments  of  EQS  must  have  the  fora  prescribed  by  MARQ  on  lines  8 and 
63: 


(X,  Y,  M,  N,  P,  E,  D,  R) 


These  eight  arguments  have  the  same  meaning  for  EQS  as  they  do  for  sub- 
routine MARQ  except  that  for  EQS,  P is  an  input  only  and  E,  D and  R 
are  outputs  only.  The  user  programs  EQS  to  obtain  — for  given  X,  Y,  M 
N and  P inputs  — ■ the  error  measure  E,  the  partial  derivative  matrix  D 
and  the  residual  vector  R. 

For  example,  the  pair  of  yaw  equations  (44)  led  to  the  subroutine 
EYAW  listed  in  Appendix  C.  The  velocity  VQ  is  passed  into  EYAW  and  the 

yaw  of  repose  is  extracted  from  EYAW  by  a labelled  CXftMON  statement 

linked  to  the  subroutine  calling  MARQ. 

The  definitions  of  the  D elements  within  EYAW  follow  immediately 
from  Equations  (44),  (48)  and  (49)  - - - with  one  exception:  I have 

ignored  the  partial  derivatives  of  the  yaw  of  repose  with  respect  to 
the  only  two  parameters  involved  therein,  and  B . Such  liberties 

may  often  be  taken  when  forming  matrix  D because  it  is  not  necessary 
that  D be  the  mathematically  correct  set  of  partial  derivatives . It 
is  possible  to  get  the  right  answer  using  a "wrong"  D array,  just  as 


we  can  get  the  right  answer  using  a "wrong"  expression  fbr  grad  c,  Eq. 
(7).  Of  course,  this  doesn't  mean  that  we  can  make  a blunder  in  coding 
D in  subroutine  EQS  or  that  we  can  write  down  any  old  approximation  that 
comes  to  mind.  A certain  amount  of  discretion  is  called  for;  if  the 
user  has  any  doubts  as  to  the  merits  of  an  approximation  in  D (and  even 
if  he  hasn't  any  doubts),  his  safest  course  is  to  avoid  such  an  approx- 
imation. 

Subroutine  MARQ  receives  matrix  D from  EQS  and  forms  the  needed 
elements  of  array  a by  summation,  in  accordance  with  Eq.  (18).  Because 
of  the  symmetry  of  a,  only  the  N«(N  + l)/2  elements  on  and  below  the 
principal  diagonal  are  computed  in  MARQ.  For  our  yaw  equations,  N is 
ten  and  thus  fifty-five  a elements  are  formed  in  MARQ.  Note  from  the 
D equations  in  the  EYAW  listing.  Appendix  C,  that  twelve  of  these  fifty- 
five  elements  should  be  identically  zero.  For  example, 

NL 

« • l [D  (I,  1)  • D (I,  3)  + D (J,  1)  • D (J,  3)] 

13  1-1 

where  J = I + NL.  Hence,  in  the  EYAW  notation, 

a - l (R1  • (-  R3)  + R3  • Rl]  - 0 
1 3 

Similarly,  the  reader  may  verify  that  except  for  sign,  only  nineteen 
of  the  forty- three  nonzero  a elements  are  distinct.  For  example, 

a • a ■ a ■ a 
62  71  84  93 

Subroutine  MARQ,  of  course,  forms  all  fifty-five  elements  by  summation; 
short  cuts  that  depend  on  the  particular  fitting  equation(s)  used  are 
sacrificed  on  the  altar  of  generality. 


XI . COMMENTS  ON  TOE  PROGRAM  THAT  CALLS  MARQ 
In  the  program  that  calls  MARQ: 

(a)  the  actual  subroutine  name  (EYAW  in  our  example)  that  will  be 
passed  as  an  argument  to  MARQ  must  be  declared  in  an  EXTERNAL 
statement. 

(b)  the  six  array  arguments  of  MARQ  must  be  declared  in  a DIMENSION 
statement: 


X (M),  Y (M),  P (N),  D (M,  N),  R (M) , EK  (N) 


Usually,  for  a given  problem,  the  value  of  N (the  number  of 
parameters)  is  known  and  fixed,  whereas  the  value  of  M (the 
number  of  data  points)  varies  from  case  to  case.  In  that 
event,  sane  number  equal  to  or  larger  than  the  largest  antici- 
pated value  of  M should  be  used  in  the  DIMENSION  statement. 

(In  our  EBL  range  set-up,  data  can  be  obtained  at  no  more  than 
54  spark  stations;  hence,  we  have  M ■ 108.) 

(c)  each  known  constant  in  the  fitting  equation  whose  value  may 
change  from  case  to  case  (such  as  the  reference  velocity  in 
our  yaw  problem)  is  assigned  a FORTRAN  name  and  passed  to  the 
EQS  subroutine  by  a labelled  COMON  statement.  Similarly,  but 
more  rarely,  any  constant  of  interest  evaluated  within  the  EQS 
subroutine  (such  as  the  yaw  of  repose)  may  be  rescued  from 
oblivion  by  linking  it  through  the  same  labelled  COMMON  with 
the  program  that  calls  MARQ. 

(d)  MARQ  will  be  called  repeatedly  (say,  in  a DO-loop)  until  some 
specified  convergence  criterion  is  satisfied  or  until  a spec- 
ified number  of  calls  have  been  made.  For  example,  we  might 
have 

DO  4 K « l.KMAX 

CALL  MARQ (EQS, X,Y,M,N,P,C, ENEW, D,R,TH,EK) 

IF(K.EQ.l)  GOTO  3 

CR  ■ 1.0  - ENEW/EOLD 
IF(CR.GE.O. .AND.CR.LT.EPS)  GOTO  5 

3 BOLD  - BNEW 

4 CONTINUE 

where  the  value  of  KMAX  (the  terminal  parameter  of  the  DO-loop) 
and  the  value  of  EPS  (where  0 < EPS  < < 1)  have  been  specified 
before  entering  the  DO-loop. 

The  process  is  assumed  to  have  converged  when  the  IF-conditions  on  CR 
above  are  satisfied,  that  is,  when  the  ratio  of  the  present  error  of 
the  fit  (ENEW)  to  the  previous  error  of  the  fit  (EOLD)  falls  between 
1-EPS  and  1.  Note  that  the  first  call  of  MARQ,  K » 1,  must  be  handled 
a little  differently  from  subsequent  calls  in  the  DO-loop  because  there 
is  no  zeroth  error  value  to  compare  with  the  error  returned  from  that 
first  call. 

Since  the  use  of  Marquardt's  X insures  that  ENEW  < EOLD,  it  may 
seem  that  CR  will  always  be  non-negative  and  hence  that  the  first  of 
the  two  IF-conditions  is  unnecessary.  The  catch  is:  whenever  we  avoid 
X (by  setting  C * - 2.0),  we  lose  our  guarantee  that  things  will  inqprove 
and  CR  can  very  well  go  negative. 

If  the  convergence  criterion  in  the  above  DO-loop  is  satisfied, 
control  is  transferred  to  statement  5.  If  (KMAX  - 1)  iterations  have 


3 


been  performed  without  satisfying  the  convergence  criterion,  then  the 
statement  following  statement  4 Is  executed  next.  What  happens  at 
these  two  locations  is,  of  course,  the  user's  affair. 

For  our  yaw  problem,  the  program  that  calls  MARQ  was  itself  written 
as  a subroutine:  the  subroutine  EY  listed  in  Appendix  D.  The  input 

arguments  are 

RD  ■ a number  identifying  both  the  range  (Aerodynamics  or 
Transonic)  and  the  round 

N ■ NL,  the  number  of  spark  stations 

NS  « an  array  of  identifying  station  numbers 

b ■ the  array:  measured  at  each  station 

A » the  array:  measured  at  «»ch  station 

Z ■ tho  array:  measured  z at  each  station,  m 

ZR  « zq,  the  reference  z,  m 

VR  ■ velocity  at  zq,  m/sec 

The  only  input/output  argument  is 

P ■ the  array  of  ten  parameters  defined  in  Eq.  (43),  with  this 

exception:  here  P(l)  and  P(6)  are  and  K^,  respectively, 

not  rij  ■ in  Kj Q (required  conversions  to  and  from  and  n2 

are  done  within  EY,  so  that  the  EY  user  need  not  be  aware 
of  the  n's).  Upon  input  to  EY,  P must  contain  the  initial 
estimates  of  the  parameters;  upon  output,  P contains  the 
final  parameter  values. 

The  output  arguments  of  EY  are: 

Q ■ the  array  of  estimated  errors  in  the  ten  parameters.  Sub- 
routine MARQ  returns  to  EY  the  estimated  errors  in  and 

n ; the  errors  in  K and  K are  then  obtained  in  EY  by 
2 10  20 
the  approximation 

E (Kj0)  - Kj0  • B (n.) 

8B  a the  array:  computed  5H  at  each  station 
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AA 


■ the  array:  computed  Gy  at  each  station 

RMS  - the  root -mean-square  error  of  the  fit  (Eq.  (20)  with  b - 

1.0) 

YREP  ■ the  yaw  of  repose,  Eq.  (42) 

IC  « a convergence  flag: 

IC  ■ 0 if  the  process  has  converged 

■ 1 if  the  process  has  failed  to  converge  in  the 
specified  number  of  iterations 

■ 2 if  the  process  has  not  been  used,  because  there 
were  too  few  (less  than  six)  spark  stations 

The  skeletal  FORTRAN  DO- loop  given  earlier  in  this  section  is 
fleshed  out  in  EY  (lines  EY  41  through  EY  57)  by  a number  of  relatively 
unimportant  statements  concerned  with  monitoring  the  progress  of  the 
convergence.  One  point  concerning  these  statements  is  of  minor  interest. 
Although  all  angles  in  array  P (that  is,  in  the  parameters  A , B , C , 

A , B and  C ) are  in  radians,  the  WRITE  statement  (line  EY  52)  prints 

2 2 2 

the  results  in  degrees;  this  is  a concession  to  the  majority  of  people 
who  "can't  picture  radians."  Someone  might  ask:  then  why  not  simply 

carry  all  angles  within  P in  degrees,  thereby  avoiding  the  need  to 
conv  :t  before  printing.  The  answer  is:  conversion  can't  be  avoided 

(unless  we  are  willing  to  print  results  in  radians).  If  the  angles  were 
carried  in  degrees,  the  conversion  that  we  eliminated  from  EY  would  crop 
up  in  EYAW,  in  taking  the  derivatives  of  sines  and  cosines.  This  would 
be  less  efficient,  since  EYAW  statements  are  encountered  more  often  in 
the  program  than  EY  statements . 

TWo  examples  of  the  print-out  furnished  by  EY  (see  Figures  3 and  4) 
are  discussed  in  the  next  section. 


XII.  A SAMPLE  CASE,  WITH  AND  WITHOUT  X 

We  define  X_ to  mean  "X  plus  the  other  features  in  Figure  2 that 
distinguish  Marquardt's  algorithm  from  the  usual  differential  corrections 
process."  If  we  apply  the  differential  corrections  process  with  and 
without  to  the  same  fitting  equation(s),  the  same  set  of  measurements 
and  the  same  starting  point  P , we  can  distinguish  three  outcomes.  We 

will  have: 

(a)  convergence  to  the  same  point,  with  or  without  X_, 
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or  (b)  divergence  without  A_,  convergence  with  X_, 

or  (c)  convergence  to  a wrong  point  without  X,  and  to  a better  (not 
necessarily  the  best)  point  with 

(A  fourth  possibility  — convergence  to  a better  point  without  X_  than 
with  it  — is  too  unpleasant  and,  I would  hope,  too  rare  to  consider 
here.) 


Note  that  when  X is  used,  divergence  is  impossible:  the  process 

can  never  blow  up.  Tf  the  terminal  parameter  of  the  DO- loop  is  large 
enough,  that  is,  if  MARQ  is  called  enough  times,  the  process  should 
converge  to  a point.  (We  need  a limit  on  the  number  of  iterations  to 
get  past  the  occasional  pathological  case.)  Of  course,  with  or  with- 
out X_,  "convergence"  in  the  above  triad  means  only  that  the  process 
has  stopped  at  some  answer;  we  may  still  be  miles  from  the  right 
answer. 

For  our  yaw  problem,  result  (a)  above  — the  least  interesting 
result  — has  occurred  by  far  the  most  often.  For  most  of  the  rounds 
fired  in  the  two  EBL  ranges,  the  motion  can  be  so  well  represented  by 
Eq.  (44)  and  initial  estimates  of  the  ten  parameters  are  so  well  deter- 
mined (by  a preliminary  subroutine  which  we  won't  discuss  here)  that 
quick  convergence  is  assured  and  Marquardt's  algorithm  is  unnecessary. 

However,  I wanted  to  complete  this  report  with  an  example  of  out- 
come (b)  or  (c),  preferably  (c)  because  I think  it  is  more  instructive. 

I could,  of  course,  have  worked  with  one  of  those  relatively  rare 
rounds  for  which  Eq.  (44)  seems  to  be  inadequate,  but  in  that  event, 
convergence  is  not  necessarily  an  advantage.  The  final  point  reached 
may  be  the  best  possible  based  on  Eq.  (44)  and  yet  be  so  worthless  that 
an  unwary  use  of  the  results  could  do  more  harm  than  no  results  at  all. 

To  obtain  outcome  (b)  or  (c)  from  any  "normal"  round,  I had  to 
resort  to  an  artifice:  I by-passed  the  subroutine  that  would  have 

given  us  good  first  estimates  of  the  ten  parameters  and  fed  mediocre 
(not  really  bad)  estimates  into  subroutine  EY  through  array  P.  In 
effect,  I said:  let's  see  what  happens  if  we  get  a little  sloppy  in 

choosing  our  starting  point.  For  round  1-11461,  the  answer  — outcome 
(c)  — is  illustrated  in  Figures  3 and  4. 

Figure  3 is  a print-out  of  the  fit  obtained  without  Marquardt's 
algorithm.  The  line  TRY  ■ 0 gives  the  initial  estimates  of  the  ten 
parameters;  these  produce  an  RMS  error  of  0.041732.  The  first  iteration 
(TRY  = 1)  improves  the  situation  but  the  second  iteration  grossly  over- 
shoots the  target.  It  is  surprising  (to  me,  at  least)  that  the  differ- 
ential corrections  process  could  recover  from  such  a flagrant  mis-fit, 
but  it  does,  converging  after  twenty  iterations  to  a point  P whose 

20 

RMS  error  is  18%  of  the  original  error. 
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Figure  4.  Yaw  fit  for  the  same  data  and  starting  point  as  in  Figure  3,  but  with 


The  point  P is  not  at  the  absolute  minimum  but  this  fact  can 
20  

hardly  be  deduced  from  a study  of  Figure  3.  Note  the  insidious  lure  of 
those  columns  of  numbers  converging  to  a wrong  answer.  If  we  had  started 
with,  say,  P as  our  initial  point,  even  our  experienced  analysts  — who 

would  suspect  any  final  result  that  differs  very  much  from  the  starting 
point  — might  be  lulled  into  acceptance.  The  error  estimates  for  the 
ten  parameters  (printed  below  the  final  iteration  values)  are  not  un- 
acceptably large;  they  help  maintain  the  illusion  that  the  parameter 
values  themselves  are  acceptable.  Finally,  the  bottom  half  of  Figure  3, 
listing  the  measured  and  computed  and  iy,  and  the  corresponding 

residuals,  seems  to  reinforce  the  impression  that  we  have  fitted  the  data 
adequately . 

And  yet  if  we  look  at  Figure  4 — the  same  round,  but  using  Marquardt's 
algorithm  — we  see  that  a much  different  fit  is  possible,  with  a much 
smaller  RMS  error  and  a much  better  set  of  parameter  error  estimates. 

(Thus  the  fact  that  the  ten  error  estimates  in  Figure  3 were  of  satis- 
factory size  proved  nothing.  If  there  had  been  no  noise,  we  would  have 
converged  there  to  some  "true  local"  minimum.  The  parameter  error 
estimates  in  Figure  3 are  measures  of  how  close  we  came  to  that  true 
local  minimum,  not  to  the  true  absolute  minimum.) 

It  is  instructive  to  see  just  how  Marquardt's  algorithm  got  us  to 
the  right  answer  in  Figure  4.  The  first  iteration  is  nearly  the  same 
as  in  Figure  3 since  X is  relatively  small  (0.001).  It  is  in  the  second 
iteration  that  Marquardt's  algorithm  had  its  first  big  opportunity  to 
star.  It  recognized  the  fact  that  a giant  over-step  was  about  to  be 
made  and  so  — since  9 is  less  than  45  degress  — it  divided  each  param- 
eter increment  by  ten.  For  example,  in  Figure  3,  the  increment  in  A 

from  the  first  to  the  second  iteration  is 

(AX.)  = - 0.18722  - 0.01878  = - 0.20600  (1/m) 

1 3 

while  in  Figure  4,  the  increment  is  about  one-tenth  as  large: 

(AX.)  = - 0.00180  - 0.01876  = - 0.02056  (1/m) 

1 4 

(For  K and  K , recall  that  it  is  the  increments  in  n and  n that  are 
10  20  12 
reduced  by  a factor  of  ten.)  A single  shrinking  of  the  step-size  was 
sufficient  in  this  instance  to  give  a smaller  RMS  error  than  obtained 
from  the  first  iteration  and  so  the  second  iteration  was  concluded. 

Thereafter  the  Marquardt  part  of  the  process  had  little  to  do. 

Usually,  it  is  not  possible  to  tell  when  the  "9  < 45°"  feature  has 
been  called  into  play  (without  inserting  monitoring  statements  within 
MARQ) . That  feature  may  or  may  not  have  been  used  whenever  the  print- 
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out  shows  a 9 less  than  4S  degrees.  We  can  say,  however,  that  X itself 
was  of  negligible  aid  in  the  convergence  of  Figure  4.  This  is  shown  by 
the  steadily  decreasing  values  of  X printed  out.  A departure  from  the 
norm 


*TRY  K “ 0,1  X 


TRY  K-l 


would  have  indicated  that  X played  a significant  role  in  completing  the 
K-th  iteration.  Such  a departure  doesn't  occur  in  Figure  4. 

Of  course,  if  we  start  out  anywhere  within  a relatively  tiny  region 
of  the  parameter  space  surrounding  the  solution  P of  Figure  4,  we  will 
converge  to  that  solution  without  X.  (This  is  what  actually  happened 
when  we  used  our  preliminary  subroutine  to  determine  the  initial  esti- 
mates.) The  point  is:  when  X^  is  used,  the  region  of  convergence  is 

expanded  considerably.  Marquardt's  algorithm  is  like  radar  on  a ship 
seeking  harbor;  on  a foggy  night,  we  need  all  the  help  we  can  get. 
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APPENDIX  A • • • SUBROUTINE  MARQ 


SUBROUTINE  MARQ(EQS.X.V»M.N.P.C.E.D.R«TH.EK» 

**** 

l 

DIMENSION  X!M)«VIM).P!Nt.D(M.NI.R(MI«EK(N)» 

MARQ 

2 

c 

1 

ALPHA! 10.101 .BETA! 10) .GAMMA (10*101 .PB(10).S(IGI 

MARQ 

3 

NOTE  — IP  NO*  CP  PARAMETERS  N.GT.IO,  REPLACE  10  ABOVE 

ANO  IN  MATINV  LIST  (LINE  MARQ  481  BY  AN  INTEGER. GE.N  . 

THE  VALUES  OF  AN.  BN  ANO  GN  BELOW  MAY  BE  CHANGED . PROVIOED 
AN  .GT.  1*.  BN  *GT.  1*  ANO  0*  «LE*  GN  *LE*  45* 

THAT 

CATA  0EC/S7. 29578/  CMIN/.5E-E6/  AN/10./  BN/10./  GN/45./ 

MARQ 

4 

PM 

■ M 

MARQ 

5 

NN 

* N 

MARQ 

6 

* GO  TO  1 tF  MARO  IS  NOT  BEING  CALLED  FOR  THE  FIRST  TIME. 

lFIC.GE.O.IGCTO  L 

MARQ 

7 

* THE  FIRST  TIME  MARQ  IS  CALLED.  EVALUATE  E*  D AND  R. 

* SET  C * INITIAL  VALUE  OF  LAMBDA  AND  RETURN. 

CALL  EQSIX.Y.MM.NN.P.E.O.R) 

MARQ 

8 

C > O.OQ1*AN*!C  ♦ 2.) 

MARQ 

9 

TH  - 0. 

MARQ 

10 

RETURN 

MARQ 

11 

* SET  CL  « INPUT  LAMBOA/AN  AND  EA  > INPUT  ERROR. 

l 

CL 

» C 

MARQ 

12 

IF 

(CL.GT.CMINI  CL  * CL/AN 

MARQ 

13 

LA 

« E 

MARQ 

14 

♦ FORM  THE  BETA  VECTOR.  EQI15I*  AND  ONLY  THE  NIN*lt/2 

* ALPHA  ELEMENTS  ON  OR  BELOU  THE  PRINCIPAL  DIAGONAL. 

EQ(I9>. 

CO 

4 J « It NN 

MARQ 

15 

6 ■ 0. 

MARQ 

16 

DO  3 K * 1(J 

MARQ 

17 

A * 0. 

MARQ 

18 

DO  2 I « 1. MM 

MARQ 

19 

IF(K.EQ.l)  B - B ♦ RIII*DII.JI 

MARQ 

20 

A « A ♦ 0( I t J >*D! I *K  1 

MARO 

21 

2 

CONTINUE 

MARQ 

22 

ALPHA! J,K)  - A 

MARQ 

?3 

3 

CONTINUE 

MARQ 

24 

BETA! Jl  - B 

MARQ 

25 

A 

CONTINUE 

MARQ 

26 
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* FORM  SCALE  FACTORS  SIJ).  REFLACK  BETA  WITH  SCALED 

* BETA*  EQ( 251.  FOAM  SCALED  ALPHAIJ.K).  EQ!25).  AND 

* STORE  ABOVE  THE  PRINCIPAL  DIAGONAL  AS  ALPHA1K.J). 

* FORM  8M  - THE  SQUARE  OF  THE  MAGNITUDE  OF  SCALED  BETA. 


C 

C 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


8M  ■ 0. 

MARQ 

27 

00  6 J » I.NN 

MARQ 

28 

S!JI  -I./SQRT1 ALPHA! J, J)) 

MARQ 

29 

BETA! J t ■ BETA  1 J )*S( J ) 

MARQ 

30 

BM  « 8M  ♦ BETA(J)*»2 

MARQ 

31 

K • J - 1 

MARQ 

32 

5 

IFIK.EQ.OIGOTO  6 

MARQ 

33 

ALPHA!  K.  J ) • ALPHA!  J.KKSI  J)*S!K) 

MARQ 

34 

K ■ K - 1 

MARQ 

35 

GOTO  5 

MARQ 

36 

6 

CONTINUE 

MARQ 

37 

* FORM  MATRIX  GAMMA.  EQ!33). 

BASED  ON  CURRENT  VALUE  OF  LAMBDA. 

7 

01  AG  « 1.  ♦ CL 

MARQ 

38 

CO  9 J ■ I.NN 

MARQ 

39 

GAMMA! J.J)  - DIAG 

MARO 

40 

K • J -1 

MARQ 

41 

8 

(F!K. EQ.OIGOTO  9 

MARQ 

42 

GAMMA! J.K 1 ■ ALPHA (K.JI 

MARQ 

43 

GAPMAIK.J)  - GAMMA! J.K) 

MARQ 

44 

K » K - l 

MARQ 

45 

GOTO  B 

MARQ 

46 

9 

CONTINUE 

MARO 

47 

* REPLACE  GAMMA  BY  ITS  INVERSE. 

CALL  MAT INVIGAMMA.NN.PBf IO.O.OOT)  MARQ  48 

* FORM  THE  COMPONENTS  OF  THE  SCALED  DELTA  P VECTOR... 

• OP  - SCALEO  (DELTA  A)  SUB  J.  SATISFYING  EQI32I. 

* FORM  THE  CANO IDA TE  POINT  PB  • P ♦ UNSCALEO  DELTA  P. 

• FORM  DOT  « OOT  PRODUCT  OF  SCALEO  BETA  AND  SCALED  DELTA  P. 

• FORM  OPM  « THE  SQUARE  OF  THE  MAGNITUDE  OF  THE  SCALED  DELTA  P. 


COT  • 0. 

MARQ 

49 

DPM  • 0. 

MARQ 

•50 

DO  11  J • I.NN 

MARQ 

51 

DP  ■ 0. 

MARQ 

52 

00  10  K - I.NN 

MARQ 

53 

CP  • DP  4 BETA !K)*GAMMAI J.K) 

MARQ 

54 

10 

CONTINUE 

MARQ 

55 

P8!J!  - P1J 1 4 OP*S!J> 

MARQ 

56 

DOT  » DOT  ♦ CP*BETA! J t 

MARQ 

57 

CPM  • DPM  ♦ OP POP 

MARQ 

58 

11 

CONTINUE 

MARQ 

59 

46 
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I 

I 

I 

I 


* FORM  ANG  - THETA*  DECREES*  E0I39). 


ANC  » 0. 

TR  ■ OOT/SQRT ( DPM4BM ) 

IF(ABS(TR)*LE*1*1  ANC  ■ OEG’ARCCOSf TR) 

* EVALUATE  E8  • ERROR  AT  THE  CANDIDATE  FOINT  PS. 

12  CALL  EQSI X* V«MM«NN*P8*EB*0*R) 

* COMPARE  ERROR  EB  AT  POINT  PB  WITH  INPUT  ERROR  EA* 

IF  4 EB.LE.EA )GOTO  19 
IFICL.EQ.O.IGOTO  19 

IFIANG.LT.  GHIGOTO  13 

* INCREASE  LAMBDA  AND  60  BACK  TO  COMPUTE  NEW  GAMMA* 

CL  - ANWCL 
GOTO  7 

* DECREASE  LENGTH  OF  OELTA  P AND  GO  BACK  TO  COMPUTE 

13  CO  14  J - i*NN 

PB(J)  « FIJI  ♦ (PB( JI-PI J|)/BN 

14  CONTINUE 
GOTO  12 

* THE  ITERATION  HAS  BEEN  COMPLETED  SATISFACTORILY. 

* UPCATE  CURRENT  POINT  P AND  ERROR  E*  COMPUTE  ERROR 

* ESTIMATES  FOR  THE  PARAMETERS* 

15  E « EB 
C • CL 
TH  > ANG 

00  16  J * 1 *KN 
PUT  - PB(J) 

EKIJI  - EB*SU)*SQRT<  GAMMA! J*J)*D( AG  ) 

16  CONTINUE 
RETURN 
ENC 


MARQ  60 
MARO  61 
MARO  62 


MARQ  63 


HARO  64 
MARQ  65 
MARO  66 


MARO  67 
MARO  68 

NEW  EB. 

MARQ  69 
MARQ  70 
MARO  71 
MARO  72 


MARO  73 
MARQ  74 
MARQ  75 
MARQ  76 
MARQ  77 
MARO  78 
MARO  79 
MARO  80 
MARQ  81 


t 


APPENDIX  R . . . SUMMIT  INK  NATINV 

OBTAINED  PROM  COMPUTER  SUPPORT  OIVISION 
ABEROEEN  RESEARCH  AND  DEVELOPMENT  CENTER 


SUBROUTINE  HATINVI A,N,C,NNAX,K,DETI 
DIMENSION  AINMAX, 1 1 »C( 1) 


IF  1 1-KK  S 3,1,1 

1 N3  - NN 

IF  IKK)  2,4,2 

2 ASSIGN  9 TO  NS 
ASSIGN  13  TO  NT 
GOTO  5 

3 N3  « KK  ♦ NN  - 1 

4 ASSIGN  1C  TO  NS 
ASSIGN  14  TO  N7 

5 GET  • 1.0 

00  15  I - l«NN 

IF  UU,m  7,6,7 

6 ViRlTE{6,17) 

OET  • 0.0 
GOTO  16 

7 T1  - 1.0/A(IfI) 

OET  - DET*A( 1 , 1 1 
All, II  - 1.0 

CO  8 J ■ 1,N3 

Afl,J)  - Al I, JI*T1 

8 CONTINUE 

GOTO  N5,  (4,101 

9 C(I)  « CIII*T1 

10  00  14  j ■ 1,NN 

IF  (I-JI  11,14,11 

11  T1  - A( J, 1 ) 

AIJ.II  - 0.0 
CO  12  L « 1,N3 

A(J,L)  • AIJ.LI  - Tl*A(I,L) 

12  CONTINUE 

GOTO  N7,  (13,14) 

13  C(J»  « C(J)  - T1*C( 1 1 

14  CONTINUE 

15  CONTINUE 

16  RETURN 

17  FORMAT  I 16H  SINGULAR  MATRIX) 

END 


*****  1 
MATINV  2 
MATINV  3 
NATINV  4 
MATINV  5 
MATINV  6 
NATINV  7 
NATINV  6 
MATINV  9 
MATINVIO 
MATINVil 
MATINV12 
MAT INVI3 
MAT INVI4 
MAT INV15 
MAT INVI6 
MAT INV17 
MATINV19 
MAT INV19 
NAT INV20 
MAT INV21 
MAT INV22 
MAT INV23 
NATINV24 
MAT INV2S 
NATINV26 
MATINV27 
NATINV28 
NAT INV29 
MAT INV30 
MATINV31 
MATINV32 
MAT INV33 
MATINV34 
MATINV35 
MATINV36 
MAT INV37 
MAT INV3B 
MAT INV39 
MATINV40 
NATINV41 
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APPENDIX  C • • • SUBROUTINE  EVAN 


SUBROUTINE  EVANIX,Y,N,N,P,E,0,AI 

pppp 

I 

(DIMENSION  X(N>»V(N),P(N)»0(P»N1,RIN1 

EVAM 

2 

COMNON/EEP/VtF 

EYAM 

3 

s ■ 0. 

EVAM 

4 

F ■ 0. 

EVAM 

5 

r - P|4)PPI9)*V*V 

EVAM 

6 

IF(T.NE.O.)  F - -9.8P|  PI41  ♦ PI9)  1/T 

EVAM 

7 

NL  - M/2 

EVAM 

8 

00  2 I "WNL 

EVAM 

9 

J • I ♦ NL 

EVAM 

10 

M • XII) 

EVAM 

II 

EL  - EXP ( PI L)  ♦ PI 2I*M  ) 

EVAM 

12 

E2  - EXPI  PI6I  ♦ PI7)*M  ) 

EYAM 

13 

AL  • PI3)  ♦ H* ( PI4)  ♦ N«PI5)  » 

EVAM 

14 

A2  • PIS)  * H* I PI9)  ♦ MPPUOll 

EVAM 

15 

RL  - EL*COSI AL I 

EYAM 

16 

R2  > £2*CGSIA2I 

EYAM 

17 

R3  • E IPS  INI  A L I 

EYAM 

LB 

R4  - E2PSINIA2) 

EVAM 

19 

RII)  - VII)  - R1  -R2  - F 

EVAM 

20 

R I J 1 « VIJ)  - R3  *R4 

EVAM 

2L 

S - S ♦ RIIIPP2  ♦ RIJIPP2 

EYAM 

22 

0(1, II  - Rl 

EYAM 

23 

01 J, L 1 - R3 

EVAM 

24 

C(  1,21  - Rl*U 

EVAM 

25 

Cl J,2 I - R3*W 

EVAM 

26 

Cl  I « 3 1 « -R3 

EVAM 

27 

D( J»  3 ! • RL 

EVAM 

28 

DU, 41  - -D( J, 2 1 

EYAM 

29 

DIJ, 4)  - DU, 2) 

EVAM 

30 

DU, SI  - -DIJ, 21PM 

EVAM 

31 

CIJ.5I  - Oil, 21PM 

EYAM 

32 

D|  1,61  - R2 

EVAM 

33 

DIJ, 61  - R4 

EVAM 

34 

Cl  1,71  - R2*M 

EVAM 

35 

DIJ, 71  » R4*M 

EYAM 

36 

ou.ai  ■ -R4 

EYAM 

37 

0IJ.8I  - R2 

EVAM 

38 

D 1 1,91  ■ -DIJ, 71 

EYAM 

39 

DIJ, 91  » 011,71 

EVAM 

40 

DII.LOI-  -DIJ, 71PM 

EVAM 

41 

D 1 J,  10  )<■  01 1 ,71PM 

EVAM 

42 

CONTINUE 

EVAM 

43 

t ■ SORT!  S/FLOAT (M- 101  I 

EVAM 

44 

RETURN 

EYAM 

45 

ENO 

EYAM 

46 

Precellm  pica  Unit 
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APPENDIX  D 


SUBROUTINE  EV 


SUBROUTINE  EY(RO,N,NS,B,A,Z,ZR,VR,P,Q,BB,AA,RHS,VREP, ICt 
EXTERNAL  EYAW 

DIMENSION  NSIN), BIN). AIN), Z<N),P<10l,X<10B),V<t08t ,01108,10), 

I R( 10S ).Qtl0t,fi8tN),AAIN) 

COPHON/EEP/VfVR 
DATA  CEG/57,29578/ 

AO  FORM ATI 1HC, 17H  TOO  FEW  STATIONS/! 

At  FORK AT  I 1H 1, 2CX, 6HR0UND  ,F8.5,20X,9IHZ*  - ,F10.4,2H  N / 

1 1H  ,54X,5HV*  ■ ,F9,3»7H  M/SEC//5H  TRY, 

2 5X,2HKl,5X,5HLAM  1,5X,2HA1,  6X,2MBl,SX,2HCl , 

3 9 X , 2HK 2 * 5X * 5HL  AM  2,9X,ZHA2,  6X,2HR2,BX,2HC2, 

A SX( AHNARQ, AX, AHNAR0,9X*  3HRMS/ 

5 1H  ,7X,2(  9X* 33HI t/Nt  (DEGt  IOEG/M)  <DEG/M**2) ,AXt , 

6 6HLAKB0A, 2X ,9H(0EG) ,9X, 9NERR0R/I 

A2  F0RKATI1H  , 14,2tF8.5,F9.5,F8.1,Fe.2,F11.6,2X) ,£9.1,F7,2,F11.6) 
A3  FORMAT! lhO.AM  ERR,2<F8.S,F9.5,F8.I,FB.2»Fll.6,2Xt ,SX* 
t 16MYAW  GF  REPOSE  - ,F8.9//t 

AA  FORMAT! ASM  •♦••WARNING****  PROCESS  FAILED  TO  CONVERGE//! 

A5  FORMAT! 5H  STA, AX, AHZIMI ,9X,?IHX1 (HI  COMP  RES1D,7X, 

1 21HXIIV)  COMP  RESI0,7X«8HDELTA  SQ,9X,AHC0KP,SX»5HRESID/I 
A6  FORMAT! IF  , I A, F 10. A, 2! AX, 3F8. A t ,AX,3F 10.6! 

NL  * N 

M « NL  ♦ NL 
DF  ■ M - 10 
IF  I0F.GT.0.I  GOTO  I 
WRITEI6, AO) 

IC  - 2 
GOTO  7 

1 1C  - 0 

V « VR 

C • -1.0 

DO  2 I * 1,  NL 

J - I ♦ NL 
XI 1)  - Z( It  - ZR 
VI  I!  » B(I) 

Y( J ) - All) 

2 CONTINUE 

Wfi  ITF(6«  Alt  RO,ZR, V 
P(l»  « ALOCI  Pit)  ) 

PI  6 ) - ALOGI  PI  6)  I 


**•* 

l 

EY 

2 

EY 

3 

EY 

A 

EY 

5 

EV 

6 

EV 

7 

EV 

e 

EV 

9 

EY 

10 

EY 

It 

EY 

t2 

EV 

13 

EY 

tA 

EY 

19 

EV 

16 

EV 

) 7 

EY 

18 

EY 

19 

EY 

20 

EY 

21 

EV 

22 

EY 

23 

EY 

24 

EY 

25 

EV 

?b 

EV 

27 

EY 

28 

EY 

29 

EY 

30 

EY 

31 

EY 

32 

EY 

33 

EY 

3A 

EY 

35 

EY 

36 

EY 

37 

EY 

38 

EY 

39 

EY 

AO 

h 


I 

I 

I 


! 


i 

i 


i 


i 


! 

i 


Printing  page  blank 
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CO  4 K > 1(26 

EV 

41 

CALL  HARQSEVAW*X( VtM( 10»P*C(E(D(R(TH(QI 

EY 

42 

ki  - expi  pm  > 

EY 

43 

A1  • CEG*PI3) 

EY 

44 

B 1 ■ CEG*PI4) 

EY 

45 

Cl  - 0EG*PI5) 

EY 

46 

*2  ■ 6XPI  PI 61  ) 

EY 

AY 

A2  ■ CEG*PI8) 

EY 

48 

B2  « CEC*P<9) 

EY 

49 

C2  • CEG*PI 10) 

EY 

50 

J - K - l 

EY 

51 

WR1TEI6«42)  JtRlfPI2)(Al(Bl»ClfR2(P(7)(A2»B2tC2(CfYH(E 

EY 

5? 

IF  IJ.EQ.C)  GOTO  3 

EY 

53 

CR  ■ 1.0  - E/EA 

EY 

54 

IF  ICR  «GE.  0.  .AND.  CR  «LT.  .0000101  GOTO  5 

EY 

55 

3 EA  - E 

EY 

56 

4 CONTINUE 

EY 

57 

IC  > 1 

EV 

58 

WRITEI6.44) 

EY 

59 

5 03  » CEG*Q( 3 ) 

EY 

60 

04  - OEG*Q( 4 1 

EY 

61 

OS  - CEG*Q( S ) 

EY 

6? 

08  « CEG*QI8) 

EY 

63 

09  • CEG*Q(9) 

EY 

64 

010-  OEG*Qt 10 ) 

EY 

65 

hi  n • ri 

EY 

66 

PI  6 1 - R 2 

EY 

67 

01 1)  - Rl*QI l) 

EY 

68 

016)  - R24QI6) 

EY 

69 

WRITE! 6, 43)  Qlll.QI 2 ) , Q3 ,04,05 ,01 6) ,0 1 7 ) ,Q8,Q9,Q10. YR 

EY 

70 

RMS  > E 

EV 

71 

YREP-  VR 

EY 

72 

WR  ITE 1 6»45) 

EY 

73 

00  6 I - I,NL 

EY 

74 

J • I ♦ NL 

EY 

75 

BBI  I ) - VII)  - RID 

EY 

76 

AMI)  - VI J ) - RtJ) 

EV 

77 

DA-  Yl  1)9*2  ♦ VI J 1**2 

EY 

78 

OB-  BBI I )**2  ♦ A A 1 1 )4*2 

EY 

79 

CC-  DA  - OB 

EY 

80 

WRITE  16. 46)  NS  II) tZ 1 f ) fY| 1 1 ,BBI 1 ) *RI 1 ) t VI J) * AAI I ) . Rl  J), 

OA.OB.OCEY 

81 

6 CONTINUE 

EV 

82 

7 RETURN 

EY 

B3 

ENC 

EY 

B4 
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LIST  OF  SYMBOLS 

A (1)  a factor  greater  than  unity  by  which  X is  to  be 

increased  if  necessary  to  insure  that  < eQ 

(2)  an  BY  input  array  argument:  the  measured  4y  values 

at  each  station 

A^,  orientation  angles  of  the  two  yaw  arms  at  zQ>  Eq.  (41); 

two  of  the  ten  yaw  parameters 

AA  an  EY  output  array  argument:  the  computed  values  at 

each  station 

AN  the  FORTRAN  name  for  A (def.  1)  in  subroutine  MARQ 

a , a , . ..a^  the  n parameters  of  the  fitting  equation  whose  values 
1 2 are  to  be  determined 

®k  *k  /°ik 

B (1)  a factor  greater  than  unity  by  which  ~tfp  is  to  be 

decreased  if  necessary  when  9 < G 

(2)  an  EY  input  array  argument:  the  measured  t'H  values 

at  each  station 


B , B 

1 l 


C , C 
1 2 


turning  rates  of  the  two  yaw  arms  at  zQ,  Eq.  (41);  two 
of  the  ten  yaw  parameters 

an  EY  output  array  argument:  the  computed  values  at 

each  station 

the  FORTRAN  name  for  B (def.  1)  in  subroutine  MARQ 
1 or  0.67449,  Eq.  (20) 

a MARQ  input/output  argument;  after  the  first  call,  C =>  X 
two  of  the  ten  yaw  parameters,  Eq.  (41) 

a MARQ  input/output  argument  and  an  EQS/EYAW  output 
argument:  the  array  of  partial  derivatives 

determinant  of  matrix  y 


55 


Jik 


LIST  OP  SYMBOLS  (continued) 

3f  (X*,  P)1 


ik 


cj 

BK 

EQS 

BY 

EYAW 

F 

f 

fH*  fV 


GN 

8 

h 

IC 


K , K 
1 2 

K , K 
10  20 


L 

M 


Dik//^ 

estimate  of  the  error  of  the  fit,  Eq.  (20);  a MARQ 
input/output  argument  and  an  EQS/EYAW  output  argument 

estimate  of  the  error  in  parameter  a^,  Eq.  (21) 

a MARQ  output  array  argument:  the  estimated  errors 

a MARQ  input  argument:  the  dummy  name  of  a subroutine 

called  by  MARQ 

the  subroutine  that  calls  MARQ  in  the  yaw  problem 
the  actual  name  of  EQS  in  the  yaw  problem 
the  criterion  function,  Eq.  (4) 
the  fitting  function,  Eq.  (1) 

the  two  fitting  functions  for  the  yaw  problem,  Eq.  (44) 

the  value  of  6 (degrees)  below  which  a new  AP  is  obtained 
by  shrinking  the  current  AP. 

the  FORTRAN  name  for  G in  subroutine  MARQ 

the  acceleration  of  gravity,  assumed  constant,  Eq.  (42) 

a dimensionless  positive  constant,  Eq.  (30-31) 

an  EY  output  argument:  a convergence  flag 

lengths  of  the  two  yaw  arms,  Eq.  (39) 

values  of  K , K at  z , Eq.  (40) 

12  0 

a local  minimum  point 

a MARQ  and  EQS/EYAW  input  argument:  the  number  of  data 

points,  m 
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MARQ 

MATINV 

m 


N 


NL 

NS 

n 

P 


RD 

RMS 


LIST  OF  SYMBOLS  (continued) 

the  subroutine  for  carrying  out  the  Marquardt  algorithm 
of  Figure  2 

a matrix  inversion  subroutine 

the  number  of  measurements:  the  number  of  data  points 

(xi*  ^ 

(1)  a MARQ  and  EQS/EYAW  input  argument:  the  numbor  of 

parameters , n 

(2)  an  EY  input  argument:  the  number  of  spark  stations, 

NL 

the  number  of  spark  stations  at  which  measurements  were 
taken,  M/2 

an  EY  input  argument:  an  array  of  station  numbers 

the  number  of  parameters 

(1)  the  parameter  set  {a^,  a^,  . . . a^};  a point  in  n- 
dimensional  parameter  space 

(2)  an  EY  and  MARQ  input/output  array  argument  and  an 

EQS/EYAW  input  array  argument:  the  paramotor  values 

a point  along  the  negative  gradiont 


the  initial  and  successive  parameter  points  in  the 
iterative  fitting  process 

the  vector  from  the  origin  of  the  parameter  space  to 
point  P 

an  EY  output  argument:  the  array  of  estimated  errors 

in  the  ten  yaw  parameters 

a MARQ  input /output  array  argument  and  an  EQS/EYAW  out- 
put array  argument:  the  residuals  of  the  fit  (observed- 

computed) 

an  EY  input  argument:  a number  identifying  both  the 

range  and  the  round 

an  EY  output  argument:  the  RMS  error  of  the  ypw  fit 
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LIST  OF  SYMBOLS  (continued) 
an  integer  > -1 

n- dimension  til  parameter  space,  in  which  the  coordinates 
of  a point  P are  the  values  of  the  n parameters 

scaled  parameter  space,  Bq.  (23),  in  which  each  co- 
ordinate has  the  dimensions  of  the  dependent  variable  y 

a MARQ  output  argument,  the  angle  6 (degrees) 

velocity  at  zQ,  m/sec 

an  BY  input  argument:  VQ;  passed  to  EYAW  (under  alias 

V)  by  labelled  COMON 

a MARQ  and  EQS/BYAW  input  argument:  the  array  of  data 

points  (x^) 

the  independent  variable  of  the  fitting  equation 

value  of  x at  which  a y measurement  was  taken, 
i ■ 1,  2,  . . . m 

the  set  of  measurements  (yi);  a MARQ  and  EQS/EYAW  input 
array  argument 

an  BY  output  argument:  SR,  the  yaw  of  repose;  obtained 

from  EYAW  (under  aliases  YR  and  P)  by  labelled  COMMON 

the  dependent  variable  of  the  fitting  equation 

the  measured  value  of  y at  x ■ x^,  i ■ 1,  2,  . . . m 

an  BY  input  argument:  the  array  of  z measurements 

an  EY  input  argument:  zq 

down-range  distance,  metres 
reference  z value,  metres 

the  curvature  matrix  H 
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Q>  Cl 


LIST  OF  SYMBOLS  (continued) 


a* 

ajk 


jk 


* 

a. 


A 


the  inverse  of  Matrix  a 
the  scaled  curvature  matrix 


ix(Sjk)? 


1 r 32  e i 

1 l**i  N, 


Eq.  (11) 


co  factor  of 
determinant  of  a 


“jk//ajj  “kk 

cofactor  of 
determinant  of 

- j (grad  e)Q,  Eq.  (13) 

* (v  ^ — 64  ■ (®i 


, t,  . . 

2 


’■V 


3e 

3\ 


, Eq.  (10) 


, Eq.  (25) 


0k//“kk 
the  modified  o matrix 


r 1 + x,  j = k 

Yjk 

^ Eq.  (33) 

L ajk,  J * k 

A 

a positive  constant  « 1 in  the 

Aa. 

i 

the  change  in  a.  in  moving  from 

J 

the  iterative  fitting  process 

Aa. 

3 

£ 

• 

cT 

< 
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I 


I 


I 


i 

I 

i 


AP 


e 


n 

2 


e 


A 


A 


M 


A 


5 


«H'  «V 


I Id 

I 30 


LIST  OF  SYMBOLS  (continued) 
the  vector  ^r0B  tIie  current  parameter 

point  to  a new  point 

F(P),  the  sum  of  the  squares  of  the  residuals  of  the 
fit  at  point  P,  Eq.  (4) 


e values  associated  with  points  PQ, 


in  (K  ),  Jin  (K  ) 

10  20 

the  angle  between  t and  AP,  that  is,  between  the  negative 
gradient  and  the  direction  actually  taken,  Eq.  (35) 

a dimensionless  positive  constant  added  to  the  diagonal 
elements  of  $ in  order  to  effect  an  interpolation  be- 
tween the  methods  of  steepest  descent  and  differential 
equations . 


the  smallest  value  of  A needed  to  satisfy  Eq.  (36) 


notation  for  the  phrase  "A  plus  the  other  features  in 
Figure  2 that  distinguish  Maxquardt's  algorithm  from 
the  usual  differential  corrections  process." 

the  complex  yaw,  ♦ i Sy,  Eq.  (39) 

the  components  of  the  complex  yaw,  Eq.  (44) 

the  yaw  of  repose,  Eq.  (42) 


orientation  angles  of  the  two  yaw  arms,  Eq.  (39) 
dimensions  of  the  bracketed  expression 


evaluation  at  the  current  set  of  parameter  values 


; 

■j 
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